Universality in the merging dynamics of parametric active contours: a study in 

MRI-based lung segmentation 
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Measurement of lung ventilation is one of the most reliable techniques of diagnosing pulmonary 
diseases. The time consuming and bias prone traditional methods using hyperpolarized H^He and 
magnetic resonance imageries have recently been improved by an automated technique based on 
multiple active contour evolution. Mapping results from an equivalent thermodynamic model, here 
we analyse the fundamental dynamics orchestrating the active contour (AC) method. We show 
that the numerical method is inherently connected to the universal scaling behavior of a classical 
nucleation-like dynamics. The favorable comparison of the exponent values with the theoretical 
model render further credentials to our claim. 

PACS: 87.10. +e, 87.68. +z, 82.60. Nh 



Ventilation analysis is an authentic way of diagnosing 
lung airway diseases. The ratio of the volume of ven- 
tilated (functional) portions of lungs to the total lung 
volume is known as lung ventilation, which is used in 
validating pulmonary drugs Q. The process involves 
two complementary magnetic resonance (MR) imaging 
modalities, the hyperpolarized helium-3 (H'^He) imagery 
and the proton (^H) imagery. Lung functionality, in- 
cluding the volume of ventilated lungs, can be obtained 
from the former modality while lung anatomic details, 
including total lung volume, are accessed through the 
latter 1] . Since manual investigation of the MR imagery 
to compute lung ventilation is extremely time consum- 
ing, an active contour (AC) or snake based automated 
method has been proposed P, Q] to compute the total 
lung volume from proton MR imagery on a 2-D slice- 
by-slice basis A snake is defined as a massless 2-D 
thin string (closed or open) that can move on the image 
domain driven by two types of forces, internal elastic 
forces and external image forces Q . Under the influence 
of these two forces an initial contour (snake) clings to 
image edges and delineates an object. A snake always 
gives continuous edges unlike any traditional edge detec- 
tor (such as the Canny method 4]), thereby eliminating 
any post-processing steps to connect the detected bro- 
ken edges. These two properties are particularly useful 
when the object outline is broken and noisy as in most 
of the MR imagery. The snake method of finding the 
object boundary relies on the initial snake placement in- 
side the image. If a small initial snake is placed inside a 
lung cavity on the MR image, while growing, the snake 
may be stopped by the associated numerical artifacts and 
may not capture the actual lung outline 0. Starting 
with a larger snake may result in missing the lung cavity 
completely. A possible solution is to start with multiple 
non-overlapping small snakes inside the lung cavity and 



evolve (grow) them until they merge with each other and 
capture the cavity outline During such a process, 
the growing snakes merge with each other into a sin- 
gle contour. This automatic merging of non-overlapping 
snakes is characterized by certain attributes: a) during 
the evolution process no two snakes overlap with each 
other, b) every snake stops evolving at the object edge 
as a single snake does during its course of evolution, and 
c) growing convex shaped snakes (e.g., circular or rectan- 
gular snakes) inside a convex object recovers the object 
boundary. Although merging of multiple snakes is ex- 
perimentally verified in a multitude of cases, a concrete 
theoretical understanding of this merging snake approach 
remains a challenge. 

The aim of this letter is to bridge this gap by showing 
that the underlying principle of multiple snake-merging is 
governed by a universal power law behavior which orig- 
inates from an inherent nucleating structure. From a 
biomedical engineering perspective, the study of the effi- 
cacy of the lung cavity delineation method is crucial for 
robust clinical application. The lung cavity segmentation 
by merging snake method involves three steps: a) initially 
small non-overlapping contours are placed inside the lung 
cavity, b) generalized gradient vector flow (GGVF) fields 
Q are computed with a Dirichlet boundary condition on 
the initial circular snakes Q], and c) all the snakes are 
evolved simultaneously and independently of each other 
with the GGVF force field as the external force for the 
snakes. This automated lung cavity segmentation is at- 
tractive for a number of reasons. While other merging 
snake algorithms, such as the one proposed by Mclner- 
ney and Terzopoulos is computationally non-trivial 
compared to the original snake evolution algorithm of 
Kass et al. Q, the merging snake algorithm by Ray et 
al. maintains the same computational simplicity of Kass 
et al. algorithm. Also, based on the position of an initial 
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FIG. 1: Evolution of multiple snakes in a 2-D slice of the 
lungs by the Ray, et al Q method. 




FIG. 2: Evolution of snakes (here circular, although one can 
use rectangular snakes as well) in a hypothetical circle image. 

snake, the rigidity parameters of the snakes can be var- 
ied, so that on one hand delicate high curvature features, 
such as costophrenic angles, can be accurately captured, 
while on the other hand, snakes can be made sufficiently 
stiff in order to avoid capturing artifacts. Fig. 1 shows 
multiple snake initialization, evolution, merging and de- 
lineation of lung boundary by the Ray et all method. 

In order to map the merging snake scenario to statisti- 
cal thermodynamical systems, we first provide the math- 
ematical background for a parametric active contour or 
snake. A snake is a curve C(s) = {p{s),q{s)) defined by 
the parameter s e [0, 1]. The snake is evolved in such a 
way that it minimizes the energy functional 01 

Es = j\^{a\C\s)f + P\C"{s)f} + E,^,[C{s)]) ds (1) 

where the first two terms give the internal energy of the 
snake (a, f3 > 0) and Ecxt represents the external energy 
added to the system. C and C" are the first and second 
derivatives of the snake with respect to s. An example of 
external force for the snake is the gradient force: E^xt — 
— |V/(x, ?/)|^, where I{x,y) is the image pixel intensity. 
In general, in the presence of an external force V(C(s, t)), 
the time dynamics of such a snake is given by jl| 

= aC"{s, t) - /3C""(s, t) + V(C(s, t)). (2) 



Note that the snake in eqn. (2) is now a function of time 
t as well. The stationary solution of eqn. (2) corresponds 
to a snake that minimizes the energy functional. Ray et 
al proposed an external force V{x, y) obtained by solving 
the following partial differential equation (PDE) applying 
Dirichlet boundary condition[l]: 

sdV/DV^V - (1 - .g(|V/|))(V - V/) = 0, (3) 

where g{a) = exp(-is:a) and f{x,y) = -i?oxt(a;, y), K 
being a tunable parameter controlling the smoothness of 
the external snake force field. Here we study the resultant 
dynamics due to the evolution of a number of such snakes, 
defined by the above system of forces. 

When multiple snakes are evolved inside the desired 
closed boundary, the growth algorithm confirms that 
they all finally merge into a single snake after a finite 
time. This merging of snakes is basically a topologi- 
cal effect and naively the phenomenology reminds one 
of " nucleation" as seen in classical first-order thermody- 
namic systems 0. In our case, we allow a finite number 
of GGVF snakes, each with a finite starting radius, to 
evolve in a 2-D plane and then numerically evaluate a 
few measurables - the nucleation time (NT), the bounded 
area (BA) after nucleation has occurred, the critical ra- 
dius ( CR ) at the time of nucleation and also the nucle- 
ation rate (NR), all as functions of snake evolution time. 
The respective quantities are defined as follows - nucle- 
ation time is the time required for all the snakes to merge 
together as a single unit, bounded area is the sum of the 
areas of the initial snakes before complete merging and 
the area under the single snake after merging, critical ra- 
dius is the radius of curvature of the merged structure 
and nucleation rate is the ratio of the number of snakes 
to the bounding area, before the merging has actually 
taken place. One should note that by complete merging 
we refer to the critical phase when all the initial snakes 
merge together for the first time. 

We perform two numerical experiments with merg- 
ing. In the first experiment, we start with a circular 
binary image of radius R containing N number of circu- 
lar snakes, each of radius r < R, randomly distributed 
inside. As described before, the initial snakes are driven 
by GGVF forces and they maintain a non-overlapping 
dynamics. In another numerical experiment, we vary the 
radii of the smaller circles and later also vary their total 
number. Both numerics are repeated with varying sizes 
of the initial domain R. The enclosed figures 3 and 4 
show the variation of the bounding area and the critical 
radius against time respectively, in non-dimensionalized 
units, in loglog plots. 

We find that all the graphs show excellent scaling in a 
broad domain, with the nonlinear zones referring to the 
saturation limits in each case. In the figures shown, 
the bounding image radius is 80 units while the cor- 
responding starting snake radius is 5 units. We have 



3 



merging mechanism of snakes we find in the AC method. 

Our starting description is that of the diffusive growth 
of two dimensional spherical droplets (circles, in 2-D) in 
a Lifshitz-Slyozov 0,13, type of continuum theory. We 
start with the thermodynamic definition of the work AW 
required to form a nucleus from an aggregate of solute 
particles 




FIG. 3: Bounding area vs time in log-log scale. 



AW = S{E - TqS + PaV) 



(4) 



where 5{E — TqS) is the increment in the free energy of 
the system and 5{PoV) is the associated external work 
done. Optimizing this total work factor, one can show 0| 
that the critical radius of a nucleating system, this being 
the rate of merger of snakes in AC method, is given by 




FIG. 4: Critical radius vs time (log-log) 
ing of the snakes. 



after complete merg- 



simulated using different number of initial snakes, indi- 
cated in the legend of the figures. The plots show power 
law variations of each of the measured quantities with 
time, i. e. y ^ and the two relevant exponents, 
that of critical radius and bounding area, have the values 
ttBA ~ 0.9 — 1.0, acR = 0.25. We have used multiple 
combinations of the parameter values but the exponents 
remain universal, that is they invariably follow a power 
law statistics. 

To analyse the results further, we begin with the hy- 
pothesis that the phenomenology of the whole process, 
that of various snakes merging together under certain 
boundary conditions in the presence of suitably defined 
external forces, is fundamentally similar to that of a do- 
main growth process as encountered in the realms of clas- 
sical nucleation. In what follows, we try to evaluate the 
time rate of evolution of a system of snakes before and 
after all the snakes have completely merged with each 
other. We start from a model of solutes trying to nucle- 
ate in a solution and then compare the results with the 
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(5) 



Here /Iq and s' are the chemical potential and molecular 
surface area of the nucleus, /3 is the surface tension and 
/i' is the chemical potential of the solute in the solution. 
As more and more solutes start nucleating i. e. snakes 
merge, the solution approaches the critical limit of su- 
persaturation and in the steady state, the critical radius 
grows as 



dm ^ dc 

dt ^dr' 



\r=R. 



(6) 



where c(i?) represents the spherically symmetric concen- 
tration distribution of snakes around a nucleus of ra- 
dius a, D being the diffusion coefficient. Following the 
treatment of Lifshitz-Slyozov 0,Q|i one can show that if 
-Rcr(O) is the critical radius at the beginning of the merg- 
ing process, then for a diffusive nature of relaxation, the 
radius follows a dynamics 



dR{t) _ Rl^jO) 



dt 



R{t) 'Rcnit) R{t) 



(7) 



We now define the dimcnsionlcss quantities x{t) = 
u{t) = -^giy and T = 2logx{t). The last one in- 
creases monotonically from to oo as the time t increases 
likewise. Combining (6) and (7) in 2D, the dynamics is 
given by 



du'^ir) 1 2 

u 



dr 



(8) 



where 7 = B.^^(Q)dx > 0. A linear stability analysis 
of the above nonlinear equation gives a fixed point at 
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We are interested in the dynamics around an 



e- neighborhood of this point 70. If 7(t) 
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(e — > as 7 ^ oo), then near the critical point {uq = 

1 /3 

(7/2) ' =3/2) the merging snakes follow a dynamics 
defined by 



du 3 2 3 2 

The time variation of the merging nuclei (snakes) is 



(9) 



27i?cR(0) 



(10) 



Thus the bounding area (~ x'^it)) of the merging snakes 
grow at the rate of t which is roughly speaking our nu- 
merical estimate (0.9-1.0) also (Fig. 3). However, the 
above analysis is only valid before complete merging of 
the snakes has occurred. In the merged phase when the 
resultant asymmetrical structure continues growing fi- 
nally to coalesce with the bounding image, the system 
dynamics is modified. To analyse the situation, we start 
from eqn.(7). The equation may alternatively be repre- 
sented as 



versa, then the snake driven by GGVF will not converge 
to the actual boundary. Although it can be shown by 
the Reed-Simon's theorem that a convex set (a cir- 
cle, say) growing within a larger convex set (a larger cir- 
cle or rectangle, say) will always merge with the outer 
boundary under the action of isotropic driving forces, to 
the best of our knowledge, no such mathematical lemma 
exists for a convex set growing in a concave set or vice 
versa. 

In summary, our achievement in the present paper has 
been to analyze the physical foundation of the GGVF 
merging technique and to provide answer to the rather 
puzzling question as to why it works so accurately. In the 
process, we have shown that the answer lies in the general 
scaling behavior of the underlying nucleation dynamics, 
defined by proper scaling laws. This clearly indicates 
that an active contour system is in the same universality 
class as the nucleation model we considered. Our results, 
in unison with probable biomedical applications, are ex- 
pected to inspire further studies in the understanding of 
lung-based diseases. 



dR{t) _ D5c{R{t)) 11 
~dr^ R(t) T mt)' 

T in the above equation is the temperature of the so- 
lution in the equivalent thermodynamic problem which 
in our case, is a measure of the average surface energy 
E, E (x T, of the system. Evidently, in 2D, E goes as 
27ri?/3. Plugging these values in the above equation, we 
see that after complete merging of the snakes has taken 
place, the effective radius of curvature of the resultant 
structure grows as i?cR(^) ^ t^^'^ (Fig- 4). Once again 
our numerical result is in exact harmony with the theory. 

Our above analysis, both numerical and analytical 
clearly suggests that below the apparently simplistic level 
of the GGVF application, the system dynamics has a 
more fundamental symmetry. This symmetry comes 
from the fact that the GGVF method lies in the same 
universality class as that of a classical nucleation model. 
This first-order critical response of the system is what 
provides the subtlety of the underlying physics. There 
is, however, a notable shortcoming with the GGVF tech- 
nique, that it does not allow us the liberty of starting 
with an initial condition at an arbitrary location. If the 
snake starts at a position in which a major portion of 
the initial snake is outside the desired boundary or vice 
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